Numerical Methods 


Erin Catto 

Blizzard Entertainment 


Sometimes the mathematical problems we are faced with in game physics are too difficult to solve exactly. Numerical 
methods provide a set of tools to get approximate solutions to these difficult problems. In particular, numerical methods 
help us get approximate solutions to the differential equations that arise in rigid body and point mass physics. 

There are many numerical methods available and it isn’t always clear which method to use in which situation. What is the 
right performance versus accuracy trade-off? Will my simulation blow up? What are the common methods used in game 
physics? 

I hope this presentation helps you better understand numerical methods and their application. 




First, a little review ... 


Differential Calculus 
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The derivative of a function is the slope 
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The derivative of a function f(x) is the slope at position x. The derivative is defined as the limit of the secant as the increment 
h goes to zero. 


Line: the derivative is constant 



function: 


derivative: 




dy 

dx 


mx + b 


= m 


The derivative of a line is the slope, which is the rise over the run. If the line is horizontal, y is constant and the derivative are 
zero. In other words, y does not change with x. 


Parabola: the derivative is a line 



We can also think of the derivative as the tangent. It is the slope at the current point on the curve. In general the slope 
changes with the x value. 


Exponential 

the derivative is another 
exponential 
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The exponential shows up a lot in physics especially when there is damping. It has the neat property that its derivative 
another exponential. 
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A wild Exponential 
Function appearded 
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You use Differentiate 
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It is not very effective 


Now a bit of Pokemon math. 


Sine wave 

: the derivative is another 
sine wave 





sine 1 

| cosine | 







function: J = sinx 
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The derivative of a sine wave is a cosine wave. 


Functions of time: the derivative is the 

rate of change 


velocity: 

dx 
x = 

dt 

acceleration: 

d 2 x 
x = . 

dt 2 


9 


In physics we usually deal with functions of time. If x is position, then x dot is velocity and x double dot is acceleration. I’ll 
be using this dot notation a lot. 


More review: differential equations 

I have the slope 
I want the original function 
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I’ve showed a few examples of you can compute a derivative given a function. What if we have the derivative and want to 
get the function back? 



A differential equation is a 
formula for the slope 


x = f(t,x) 

(first order differential equation) 
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This is the standard form of a first order differential equation in time. 


Solving the differential equation 
means getting the original function 


x = f(t,x ) 
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Solving a differential equation means getting the original function back. If our function f depends only on t, we only need to 
perform integration. However, if f depends on x, then solution becomes more challenging. 



The simplest differential equation 


x = 0 


solution: x = X 0 

(constant) 
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So lets start solving some differential equations! 


Here’s is the simplest differential equation I could find and its solution. 


This is just a line. 


Another differential equation 



solution: x = at + X Q 


14 


What if x is on the right? 


x — ax 


solution: 


X = x 0 e 


at 


d 


proof: — (x 0 e at ) = ax 0 e al = a(x 0 e a '^ = ax 
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The exponential is needed to solve this differential equation. 
Remember that the derivative of an exponential is another exponential. 


Newton’s 2nd Law 


ma = F 


or 

nix = F 


(second order differential equation) 
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Newton’s law is a differential equation for the acceleration. Because it involves the second derivative, it is a second order 
differential equation. 




Solution 


• • 

y = -g 

y = v 0 -gt 

1 2 

y = y 0 +v 0 t--gt 
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Solving for a constant force just requires the integration (anti-derivative). I’m also accounting for the initial position yO and 
velocity vO. 




Parabolic motion 



This is why projectiles follow a parabolic path. 



Mass Spring 


x 



ground 


mx + kx = 0 
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the harder it pulls. 


The mass spring system has a force that depends on position. The more the spring is stretched, 
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Angular frequency 


x + co 2 x = 0 



radians per second 
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I can divide through by the mass to go from two parameters (m,k) to one parameter (omega). 

Omega is the angular frequency of the mass-spring system. Notice that the angular frequency increases with the spring 
stiffness and decreases with the mass. 




Solution 


x = — sintf# + x 0 cos&tf 
co 
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The solution is a sine wave with angular frequency omega. We need to combine a sine wave and cosine wave to account for 
initial conditions. 




Position and velocity versus time 



time 



time 
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v 0 =0 
(0 = 1 


These graphs show position and velocity versus time. Notice that position and velocity are 90 degrees out of phase. So the 
position is a cosine wave and the velocity is a sine wave. 


Phase Portrait: velocity versus postion 



x 0 = l 
v 0 = 0 
CO = 1 


A nice way to visualize the oscillation is to plot velocity versus position. You can view this as a parametric curve with time as 
the curve parameter. This is drawn here as a unit circle, but keep in mind that position and velocity have different units. The 
arrows show the direction of increasing time. The position and velocity start at the green circle and progressing clockwise 
around the circle. And the state continues around the circle forever. The angular frequency controls how fast the state goes 
around the circle as time moves forward. 

This graph is sometimes called the phase plot or phase portrait. This plot shows all possible states (phases) of the system 
as time goes forward. 



This animation shows position and velocity versus time along with the phase portrait. 





System of differential equations 



So far I have showed some simple physical systems, such as a projectile and a mass-spring system. 

The simulations in modern games are more complex. Often there are multiple rigid bodies colliding with friction. This means 
that there are multiple mutually dependent differential equations. We have little hope of finding the exact solutions for these 
complex systems of differential equations. 


Therefore we need to turn to numerical methods. 



Numerical solution 


• Sacrifice accuracy 

• Time stepping 

• Generic simulation 
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We can simulate systems of differential equations using numerical methods. This means we are abandoning exact answers. 
But we are making games, so this is okay. You will see that numerical solution works well with typical game loops. 
Numerical solution also creates a nice abstraction that allows us to solve differential equations in a generic manner. These 
methods can handle any number of rigid bodies and any number of forces (as memory and processor power allow). 



Generic differential equation 


f ^ 

x 


u 


x = < 


> 


x„ 



first order system of differential equations 
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Numerical methods like to have differential equations in a standard form. We can have multiple differential equations 
assembled into an array. That way we can solve them all at once. This also allows us to account for objects that are 
interacting and have mutually dependent solutions, such as stacking of rigid bodies. 


Handling second order differential 

equations 


have: mx = F 

want: X = f(x) 
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I will be dealing with Newton’s law and I need a way to get this into the first order form needed by numerical methods. 


Split into two differential equations 


. F 

V = — 

m 
x = v 



v 

x 


F_ 

m 

v 


X 
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The idea is to convert a single second order differential equation into two first order differential equations. This is 
accomplished by introducing the velocity as an independent state. Then I get a first order DE for the velocity and a first 
order DE for the position. 



Finite Difference 
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The basic approach of numerical solution is to use the finite difference. This allows me to convert a differential equation into 
an algebraic equation. There are many kinds of finite difference formulas. This one is the forward difference. It approximates 
the derivative as the next value x2 minus the current value xl over the time step h. 

You can see by the green line that this really is an approximation of the slope at xl . The approximation gets better by 
decreasing the time step h. However, taking smaller time steps is more expensive because you need to take more of them. 


Backwards difference 



We can also flip things around and approximate the derivative based on the previous value. This is called the backwards 
difference. 


Numerical solution via explicit Euler 


X = 





— L = f{ t l’ x l) 


x 2 =Xj +hf[t 1 ,x l ) 


explicit Euler 
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Lets take our standard first order differential equation and apply the forward difference formula to remove the derivative. 
Notice that this converts the differential equation into an algebraic equation. This leads to the explicit Euler method: we can 
simply plug in the current value for xl and tl to compute x2. 



I can apply explicit Euler to Newton’s law. Velocity is extrapolated based on the current force and position is extrapolated 
based on the current velocity. 




Pseudocode 


struct State 

{ 

float t , x, v; 

} 

State ExplicitEuler ( float h. State statel) 

{ 

float force = ComputeForce (statel) ; 

State state2; 

state2.t = statel. t + h; 

state2.v = statel. v + h * force / mass; 

state2.x = statel. x + h * statel. v; 

return state2; 

} 
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Explicit Euler for mass-spring system 


r\ 

v 2 = v, - hco x, 
x 2 = x l + hv { 
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Here is explicit Euler for the mass-spring system. So here omega is the angular frequency in radians per second. Notice how 
explicit Euler gives us a way to simulate an oscillation without using trig functions. 




Explicit Euler phase portrait 



Here is a Matlab simulation of the mass-spring system using explicit Euler. As before it is a plot of velocity versus position. 
The exact solution is a circle. But explicit Euler is diverging and the amplitude is growing unbounded as time progresses. 


Explicit Euler is extrapolation 
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Explicit Euler is essentially extrapolation. Intuitively extrapolation is bad. If the function is rapidly changing, the extrapolation 
has poor accuracy and can lead to a numerical instability. 


Stability of explicit Euler 


V 2 

*2 


+ h 


-co 2 x l 

V 1 
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The Matlab simulation showed that explicit Euler is unstable. We can also understand the instability by looking at the explicit 
Euler equations in detail. For the mass-spring system, the time stepping process is equivalent to multiplying the previous 
state vector by a matrix A. 


Stability is determined by growth 


x 2 = Ax, 

x 3 = A 2 x, 


X 


n 



stable if: 


A <1 
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As time marches forward, I am just multiply the previous state by the same matrix over and over. The magnitude of this 
matrix must be less than one, otherwise the solution will blow up. 


Use eigenvalues to determine 
the magnitude of the matrix 


A =A 


max 


Au = Au 


eigenvector 


eigenvalue 
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The magnitude of a matrix can be represent by its eigenvalue lambda, which is a scalar. Please see Wikipedia for more 
information on eigenvalues. 


Use the determinant to compute the 

eigenvalue 


(Al-A)u = 0 
det(AI- A) = 0 
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I can compute the eigenvalue of a matrix using the determinant. 




Eigenvalue for explicit Euler 


det 

f 

A 

1 

0 


1 - hco 2 



0 

1 


h 1 


A 2 - 2A + 1 + h 2 co 2 = 0 characteristic polynomial 
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The determinant produces a polynomial in lambda. This is called the characteristic polynomial. You can compute the roots 
of the polynomial to determine the values of lambda. In this case there are two values of lambda, but I just care about the 
value with the largest magnitude. 


Eigenvalue for explicit Euler 


X - 1 ± ihco 


X\ = yjl + h 2 CQ 2 



(unstable if h and omega > 0) 
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The eigenvalue is a complex number, but its magnitude is greater than one if the spring stiffness is positive. This shows that 
explicit Euler is unstable. You can get explicit Euler to be stable by adding some damping, but as you will see, there are 
better alternatives. 



Numerical solution via implicit Euler 


x = f(t,x) 




backwards difference 


r 
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Now lets use the backwards difference instead of the forward difference. After shifting the indices, this leads to an implicit 
function for x2 and the method is called implicit Euler. To actually compute x2 we have to solve an algebraic equation. It may 
be nonlinear, so we may have to use numerical root finding, such as the Newton-Raphson method. 



Implicit Euler for mass-spring system 


2 

v 2 = Vj — hco x 2 
x 2 = x { + hv 2 
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Here is implicit Euler for the mass-spring system. Notice that x2 and v2 appear on both sides of the formulas. This means I 
have to solve for x2 and v2 to advance time. 

I can solve these equations easily since they are linear and only have two variables. Larger simulations and non-linearities 
can make implicit Euler much more challenging to implement. 




Implicit Euler phase portrait 



The phase portrait shows that implicit Euler is stable. Also, implicit Euler is damped. So even though the mass should 
oscillate forever, implicit Euler damps the motion and brings the mass to rest. I call this numerical damping. 

If you compare this figure to explicit Euler, you should see that they look they look similar, except that explicit Euler is 
growing. This is because implicit Euler gives the same result as running explicit Euler in reverse. 


Implicit Euler eigenvalue 




l 

\l I +lror 



(always stable) 
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I performed a similar computation for implicit Euler and the result shows that implicit Euler is stable for any positive value of 
omega. 



Semi-implicit Euler: cheap and stable 


Explicit 


Semi-implicit 


F, 


v 2 = v i +h — 


m 

x 2 = x l + hv l 


v 2 = v l +h — 
m 

x 2 =x ] +h@ 
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Explicit Euler is efficient but has poor stability. Implicit Euler has great stability but can be expensive. There is a 3rd method 
that strikes a balance between these two methods, giving good stability and performance. 

If I use the current velocity in the position update, then explicit Euler becomes semi-implicit Euler. This small difference 
makes a huge difference in stability. Yet the computational cost is just as good as explicit Euler. 

Semi-implicit Euler is often called symplectic Euler because it has some good energy conservation properties. You’ll see that 
shortly. 


Semi-implicit Euler eigenvalue 


A 2 +(h 2 co 2 -2)A + 1 = 0 characteristic polynomial 


|A| = 1 if hco< 2 


conditionally stable 
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An eigenvalue analysis shows that semi-implicit Euler is conditionally stable. This is much better than explicit Euler and 
practice we can push semi-implicit Euler quite hard before it goes unstable. 



Semi-implicit Euler is conditionally stable 



At omega = 1 , semi-implicit Euler is stable and accurate. Additionally it conserves energy on average. 


ncreasing the frequency 



position 


Increasing the angular frequency reveals some deviation from the exact solution. But it is still stable and still conserves 
energy on average. 


co- 4 






ncreasing the frequency 



co = 10 


/j = 0 . 1 


hco = l 


Half the stability limit, things get strange, but it is still stable. 



(0 = 10 






ncreasing the frequency 



o 

posifop 


30 


20 


-30 


£0 = 16 


h = 0A 


This get really wacky at larger frequencies. Yet it didn’t blow up yet. 


ncreasing the frequency 



This shows that eventually semi-implicit Euler will blow up. 

Rule of thumb: aim for at least 4 time steps per spring oscillation. 


Your first choice should be 
semi-implicit Euler 

Runge Kutta 

Implicit Euler 




Semi-implicit Euler 


Verlet 

Mid-point 


Explicit Euler 
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So far I have shown 3 numerical methods for differential equations: explicit Euler, implicit Euler, and semi-implicit Euler. 

There are many more algorithms out there, such as Runge-Kutta, mid-point, predictor-corrector, and more. For games, your 
first choice should be semi-implicit Euler. It is fast, reasonably stable, and easy to implement. 

The increased accuracy of Runge-Kutta, such as RK4, may seem appealing. However, any simulation with rigid body 
collisions is already wildly inaccurate. What matters more is stability and performance. In this department semi-implicit Euler 
delivers. 


Application to games 


Gyroscopic Motion 
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Now I would like to show an example of applying numerical methods to solve a difficult problem in rigid body dynamics, 
particular we are going to look at gyroscopic motion of rigid bodies. 



Rigid body rotation 


Icb + coxIco = T 

I : inertia tensor 
co : angular velocity 
cox Icq : gyroscopic torque 
T : applied torque 


co 



L = Ico 

angular momentum 
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This equation is known as Euler’s equations of motion. You can view it as Newton’s law for rotation. The inertia tensor is 
represented by a 3-by-3 matrix. The angular velocity and the applied torque are both 3D column vectors. Note that now 
omega represents angular velocity, not angular frequency! Curiously, they both have the same units: radians per second. 

This equation looks like ma = F, except there is a middle term that depends on angular velocity. This is the gyroscopic 
torque. It accounts for non-spherical shapes that may wobble as they spin freely in the air. Specifically, it arises because 
angular velocity and angular momentum may point in different directions. 




The gyroscopic torque is due 
to a non-uniform inertia tensor 



The inertia tensor for a box with non-equal sides is non-uniform. So apply the inertia tensor to the angular velocity results in 
an angular momentum that points in a different direction than the angular velocity. So the cross product between the angular 
velocity and the angular momentum is non-zero. 


Semi-implicit Euler 


co 2 = co x + hl x 1 (T x - co jX I l co l ) 
h 

Qi ~ 2 
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This is semi-implicit Euler applied to the rotational equation. I’m using a special formula to update the quaternion. See the 
references for details. 




Capsule spin test 
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Here is a simulation of a spinning capsule using semi-implicit Euler. This video was created using the Blizzard physics 
engine called Domino. 



Why is this unstable? 


co 2 = co l + hl x 1 (T x 


co x x I i co l 


quadratic in omega 
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Here is the velocity update from semi-implicit Euler. This blows up because semi-implicit Euler extrapolates velocity. Since 
the gyroscopic torque is quadratic in the angular velocity, the extrapolation easily diverges. 


Typical trickery 


Drop gyroscopic term: Ico = T 


Angular momentum as state: 


L = I co 
L = T 
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There are a couple tricks used to avoid the gyroscopic term. Many physics engines simply drop the gyroscopic term. This is 
what I have done for many years. 

Some folks like to avoid the gyroscopic term by using the angular momentum L rather than the angular velocity. This makes 
the gyroscopic term disappear from the rigid body equations of motion. 


The problem with angular momentum 


Z^2 — L\ + hT 

®2 = | A | ^2 

h 

0.2 0 \~^~ “ ^ 2^1 


old inertia tensor : incorrect angular velocity 
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I can propagate the angular momentum with good stability using semi-implicit Euler. I need the updated angular velocity to 
compute the new orientation (quaternion). To compute the new angular velocity I need to apply the inverse inertia tensor to 
the angular momentum. In general the inertia tensor depends on orientation, so I am forced to use the old inertia tensor. This 
loses the correct rotational behavior. 


Lack of gyroscopic term in Diablo 3 
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In the past I’ve run into real game scenarios where dropping the gyroscopic term caused some odd behavior. In this case 
this creature’s staff is a slender box that spins a bit too much. 




Dropping the gyroscopic term 
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The capsule test is now stable. But it is too stable! The capsule spins forever. 



. .. adding rolling resistance 
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Adding rolling resistance results in a simulation that is still stable, but not plausible. 



Superposition Principle 

stage 1 explicit: A co [ = hI~ 1 T ] 

stage 2 implicit: A co 2 = -hl ~ 1 ( co 2 x I 2 co 2 ) 

cq 2 =cq { + A co l + A (o 2 


I can include the gyroscopic torque by using superposition principle which says that the effect of multiple separate torques 
is additive. This lets me use implicit Euler on the gyroscopic torque while using explicit Euler on all the torques. 



Solving the implicit equation 


I 2 ( co 2 — C 0 \ ) + hco 2 x I 2 co 2 = 0 


f(to 2 ) = 0 


72 


For implicit solution of the gyroscopic term, I need to solve for omega. However, this is a nonlinear equation in 3 variables 
(xyz rotation). So I solve this using Newton-Raphson iteration. Therefore I put the problem into a standard form of some 
function of omega equals zero. 




Problem: inertia tensor 


I don’t have / 2 


Solution: work in local coordinates 


The inertia tensor is constant in local coordinates 
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The inertia tensor in world coordinates depends on the orientation of the rigid body. So we cannot compute 12. By working 
in local coordinates, I can side-step this problem. 


Newton-Raphson solver 


m = o 



single variable 


f (x) = 0 

J( x „ + i- x J = f „ 

j_ df_ 

a.x, ■" dx n 
multiple variables 
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Here is a summary of Newton-Raphson. Typically it is described as a method for solving a single scalar equation. However, 
it can be extended to solve multiple simultaneous equations. To do this I need to compute the Jacobian of the function. 





Normally the inertia tensor is a function of orientation. However, I can build the implicit function for omega in local 
coordinates (b for body coordinates). The inertia tensor is constant in local coordinates so it does not need to be 
differentiated when I compute the Jacobian. Here the skew function takes a vector and builds the corresponding skew 
symmetric matrix. This vastly simplifies the solution. 




Gyrocopic torque solver 


Vec3 SolveGyroscopic (Quat q. Matrix lb, Vec3 omegal, float h) 

{ 

// Convert to body coordinates 

Vec3 omegab = InverseRotate (q, omegal); 

/ / Residual vector 

Vec3 f = h * Cross (omegab, Multiply33 ( lb, omegab)); 

// Jacobian 

Matrix J = lb + h * (Multiply33 (Skew (omegab) , lb) - Skew (Multiply33 ( lb, omegab))); 

// Single Newton-Raphson update 
omegab = omegab - Solve33(J, f ) ; 

// Back to world coordinates 
Vec3 omega2 = Rotate (q, omegab) ; 
return omega2; 
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Here is the implicit gyroscopic force solver. It uses body coordinates because I do not have the new world inertia tensor yet 



Capsule implicit solution 
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Here is the result of using implicit solution of the gyroscopic torque. I only used one Newton-Raphson iteration. 



With rolling resistance 
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Adding rolling resistance makes the solution even better. 



The Dzhanibekov Effect 



Here is a video of the Dzhanibekov effect shown on the ISS. Similar effects can be seen in tennis rackets, white board 
erasers, and mobile phones (when they slip out of our hands, that is why they are so hard to catch). 


Explicit gyroscopic integration 
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I can the experiment using explicit Euler. Unfortunately the kinetic energy grows unbounded. 



Dropping gyroscopic term 



81 


Here I dropped the gyroscopic term and you can see there is no tumbling motion. 




In this video, I have included the implicit gyroscopic solver using one Newton-Raphson iteration. The tumbling motion is 
reproduced. You can see that the kinetic energy dissipates over time. 

This sort of tumbling motion can make destruction and explosions look better in games, especially when the are long skinny 
objects. 



Conclusion: 

numerical 

qualities of 
methods 

• Convergence 


• Accuracy 


• Stability 
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There are few qualities of numerical methods that are used to find the best match for what our games need. 

Convergence: as the time step gets smaller, the numerical result approaches the exact solution. For example, a numerical 
method that always returns 42 is not convergent for a harmonic oscillator. 

Accuracy: how well the solution matches the solution for a given time step. In games we can set the accuracy bar low. We 
aren’t building satellites. 

Stability: the numerical result grows unbounded even though the exact solution is unbounded 
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Appendix 


Here is stuff that didn’t fit into a one hour presentation. 



Mass-Spring-Damper 


x 

spring k 


mass m 


ground damper c 
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mx + cx + kx = 0 


It is also interesting to add damping to see the effect on the phase portrait. 


Rewrite 


mx + cx + kx = 0 
x + 2,^cox + co^x = 0 



angular frequency 



damping ratio 
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We can reduced from 3 parameters (m,k,c) to 2 parameters (omega, zeta) by dividing through by the mass. The new 
parameter zeta is called the damping ratio. 

If zeta is zero, there is no damping. Increasing zeta increases damping. Once zeta reaches 1 , the damping is so strong that 
all oscillation is eliminated. 



Solution 


x = e t ’ m (Cj sin co d t + C 2 coso>/) 


c o d = co^l l-^ 2 damped frequency 
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Here is the oscillatory solution to the mass-spring-damper. The constants Cl and C2 depend on the initial conditions. 


When zeta is larger than one, then this solution no longer works. In that case the solution is just a decaying exponential. 


Phase portrait 
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Here is the phase portrait with a small amount of damping. As the time moves forward the position and velocity both decay 
to zero until the mass comes to rest. 



Phase portrait animation for mass-spring-damper system 





Similarity of explicit and implicit Euler 


explicit 


implicit 
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The phase portrait shows that implicit Euler is stable. Also, implicit Euler is damped. So even though the mass should 
oscillate forever, implicit Euler damps the motion and brings the mass to rest. I call this numerical damping. 

If you compare this figure to explicit Euler, you should see that they look they look similar, except that explicit Euler is 
growing. This is because implicit Euler gives the same result as running explicit Euler in reverse. 





Duality of explicit and implicit Euler 


=X i +hf i 

explicit 
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x l+l =x,+hf i+l 

implicit 


This shows that we can run time backwards on explicit Euler and get the same formula that is used for implicit Euler. 



Taylor Series 


1 2 1 3 

x 2 = x { + hx x + -Hj + — h 3cj + . . . 

2 6 
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Finite difference formulas come from the Taylor Series, which states that any function can be extrapolated by using its 
current value and derivatives. The more terms we use, the more accurate the extrapolation. If we include the first derivative 
of x, we have first order accuracy. Including the second derivative gives second order accuracy, and so on. Notice that when 
h is small, high order derivatives contribute decreasing amounts to the approximation. 




Taylor Series: sine wave 


^ = sin(r) 



increasing 

accuracy 


Here is the Taylor Series approximation of the sine function with increasing accuracy as the number of terms increases. Note 
that the approximations are polynomials in time. 

The time could be just our time step h because we can start the Taylor Series extrapolation at any point of time. 


Taylor Series: sine wave 

O Ntef ■ 0 
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Here is how the Taylor Series can approximate a sine wave with increasing order. 



First Order Approximation 



~ x, + hx x 



x i = 


X 2 -Xj 


h 


truncated Taylor Series forward difference 


96 


The first order approximation of the Taylor Series only keeps the first derivative. This is what I used already for the forward 
difference. 


Second Order 




2 


~ x, + hx, + — h 

1 1 2 



problem 
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Keeping the second derivative results in a second order approximation. I can use this to derive a second order finite 
difference approximation. However, it may not be easy to compute the second derivative of x when x defined by a 
differential equation. Fortunately, there is a trick for eliminating the second derivative. 


Central Difference 


TS1 

TS2 


x 0 = x, - hx x + — h 2 x x 


1 


x 2 = Xi + hx x + — h 2 x x 


TS2 - TS1 


x 2 - x 0 = 2 hx 


i 




Here I have two second order Taylor Series approximations. By subtracting the first equation from the second equation, I 
can eliminate the second derivative. This results in the second order finite difference formula called the central difference. 
This shows the power of the Taylor Series. You can combine more Taylor Series instances to generate ever more accurate 
finite difference formulas. Although this is exciting, it is not really necessary for games. We can get by on first order accuracy 
most of the time. 


Central Difference 



Here’s how the central difference looks graphically. It really does look like a better approximation of the slope at tl . 


